با استفاده از داده های زلزله ها در ایران و جهان به سوالات زیر پاسخ دهید.
توابع، دادهها و کتابخانههای زیر در طول تمرین مورد استفادهقرار میگیرند.
library(readr)
library(dplyr)
library(stringr)
library(ggmap)
library(highcharter)
library(plotly)
library(ggplot2)
library(rworldmap)
library(countrycode)
library(cowplot)
library(sp)
library(gganimate)
library(lubridate)
earthquakes = read_rds("data/historical_web_data_26112015.rds")
# convert coordinate to country
coords2country = function(points)
{
countriesSP = getMap(resolution='high')
pointsSP = SpatialPoints(points, proj4string=CRS(proj4string(countriesSP)))
indices = over(pointsSP, countriesSP)
indices$ISO3
}۱. با استفاده از داده های historical_web_data_26112015.rds و استفاده از نمودار پراکنش سه بعدی بسته plotly نمودار طول، عرض و عمق زلزله ها را رسم نمایید. علاوه بر آن بزرگی هر نقطه را برابر بزرگی زمین لرزه قرار دهید.
plot_ly(earthquakes,
x = ~Longitude, y = ~Latitude, z = ~Depth, size = ~Magnitude,
marker = list(color = ~Magnitude,
symbol = 'circle',
sizemode = 'diameter',
colorscale = c('#FFE8A8', '#884848'),
showscale = TRUE),
sizes = c(0.8, 8.8),
text = ~paste('Province: ', Province,
'<br>City: ', City,
'<br>Earthquake Magnitude: ', Magnitude)) %>%
add_markers() %>%
layout(scene = list(title = 'Iran Earthquakes',
xaxis = list(title = 'Longitude', range = c(40, 80)),
yaxis = list(title = 'Latitude', range = c(20, 40)),
zaxis = list(title = 'Depth', range = c(0, 180))))۲. پویانمایی سونامی های تاریخی را بر حسب شدت بر روی نقشه زمین رسم نمایید.(از داده زلزله های بزرگ استفاده نمایید.)
disaster = read_delim("data/disaster.txt",
"\t",
trim_ws = TRUE) %>%
select(lon = LONGITUDE,
lat = LATITUDE,
mag = INTENSITY,
dep = FOCAL_DEPTH,
country = COUNTRY,
year = YEAR,
t_death = TOTAL_DEATHS) %>%
na.omit()eq_an = ggplot() +
borders("world", colour="black", fill="gray") +
geom_point(aes(x = disaster$lon,
y = disaster$lat,
size = disaster$mag,
frame = disaster$mag),
color="blue") +
ggtitle("Earthquakes In The World") +
xlab("Longitute") +
ylab("Latitude")
gganimate(eq_an, "eq_an.gif") ***
۳. نمودار چگالی دو بعدی زلزله های تاریخی ایران را رسم کنید.( از داده iran_earthquake.rds و لایه stat_density_2d استفاده نمایید).
eq_ir = read_rds("data/iran_earthquake.rds") %>%
filter(Long < 70)
ggplot(eq_ir, aes(x=Long, y=Lat)) +
geom_point() +
stat_density_2d(aes(fill = ..level.., colour = ..level..), geom = "polygon") +
ggtitle("IR Earthquakes Density") +
xlab("Longitude") +
ylab("Latitude")۴. احتمال اینکه در ایران در پنج سال آینده زلزله به بزرگی هفت ریشتر رخ دهد را محاسبه کنید. (از احتمال شرطی استفاده کنید.)
برای این سوال احتمال اینکه در پنجسال آینده زلزله رخدهد به شرط آنکه بیشتر از ۷ ریشتر باشد را محاسبه میکنیم. برای اینکار احتمال اینکه در پنجسال آینده یک زلزله رخدهد را بهشرط آنکه آن زلزله بزرگتر از ۷ ریشتر باشد را بایستی محاسبه کنیم. در پنجسال آینده باتوجه به دادهها با احتمال ۱۰۰٪ زلزله رخخواهدداد پس بایستی فقط احتمال اینکه در پنج سال آینده زلزلهی به بزرگی بیشتر از ۷ ریشتر رخدهد را محاسبه کنیم.
eq_ir_7 = disaster %>%
filter(mag >= 7, country == 'IRAN') %>%
arrange(-year) %>%
mutate(dis = lag(year),
dis = dis - year,
dis = ifelse(is.na(dis),2018-year,dis))
eq_ir_7_under5 = eq_ir_7 %>%
mutate(uder5 = ifelse(dis <= 5, 1, 0)) %>%
group_by(uder5) %>%
summarise(count = n())
# +3 for continues < 5 years
prob = eq_ir_7_under5 %>%
summarise(prob =((7 + 3)/sum(count)))
knitr::kable(prob, 'html')| prob |
|---|
| 0.625 |
۵. بر اساس داده های زلزله های بزرگ ابتدا تعداد و متوسط کشته زلزله ها را بر حسب کشور استخراج نمایید. سپس نمودار گرمایی تعداد کشته ها را بر روی کره زمین رسم نمایید.(مانند مثال زیر!)
disaster_country = disaster %>%
group_by(country) %>%
summarise(death_count = sum(as.integer(t_death)),death_mean = mean(as.integer(t_death))) %>%
arrange(-death_count) %>%
mutate(country_code = countrycode(country,
"country.name",
"iso3c"))
knitr::kable(disaster_country %>% select(-country_code))| country | death_count | death_mean |
|---|---|---|
| CHINA | 609347 | 18465.06061 |
| IRAN | 298903 | 15731.73684 |
| ARMENIA | 187000 | 62333.33333 |
| ITALY | 149711 | 7129.09524 |
| PAKISTAN | 146342 | 24390.33333 |
| TURKMENISTAN | 110001 | 55000.50000 |
| PERU | 83055 | 2373.00000 |
| AZERBAIJAN | 80087 | 26695.66667 |
| ECUADOR | 71202 | 17800.50000 |
| CHILE | 65469 | 2618.76000 |
| TURKEY | 59062 | 4218.71429 |
| INDIA | 55233 | 6904.12500 |
| VENEZUELA | 27808 | 3972.57143 |
| GUATEMALA | 23054 | 4610.80000 |
| TAJIKISTAN | 15529 | 5176.33333 |
| MOROCCO | 13728 | 6864.00000 |
| RUSSIA | 12025 | 1717.85714 |
| NEPAL | 11771 | 3923.66667 |
| PHILIPPINES | 10883 | 473.17391 |
| MEXICO | 10349 | 862.41667 |
| ALGERIA | 7308 | 1461.60000 |
| JAPAN | 6125 | 322.36842 |
| GUADELOUPE | 5000 | 5000.00000 |
| UZBEKISTAN | 4880 | 4880.00000 |
| GREECE | 3133 | 120.50000 |
| YEMEN | 2800 | 2800.00000 |
| TAIWAN | 2325 | 332.14286 |
| INDONESIA | 2202 | 100.09091 |
| COLOMBIA | 2158 | 269.75000 |
| EL SALVADOR | 1328 | 265.60000 |
| MACEDONIA | 1083 | 361.00000 |
| AFGHANISTAN | 934 | 155.66667 |
| USA | 655 | 27.29167 |
| ARGENTINA | 516 | 103.20000 |
| KAZAKHSTAN | 451 | 225.50000 |
| GUINEA | 443 | 443.00000 |
| GEORGIA | 280 | 93.33333 |
| NEW ZEALAND | 258 | 86.00000 |
| ALBANIA | 243 | 121.50000 |
| PAPUA NEW GUINEA | 180 | 30.00000 |
| USA TERRITORY | 168 | 84.00000 |
| MONTENEGRO | 131 | 131.00000 |
| BULGARIA | 130 | 65.00000 |
| PORTUGAL | 113 | 113.00000 |
| COSTA RICA | 93 | 31.00000 |
| KYRGYZSTAN | 79 | 39.50000 |
| ETHIOPIA | 70 | 35.00000 |
| AZORES (PORTUGAL) | 69 | 69.00000 |
| SOLOMON ISLANDS | 62 | 20.66667 |
| MONGOLIA | 30 | 30.00000 |
| CANADA | 28 | 28.00000 |
| GHANA | 22 | 22.00000 |
| SLOVENIA | 22 | 11.00000 |
| IRAQ | 20 | 20.00000 |
| ROMANIA | 18 | 6.00000 |
| AUSTRALIA | 12 | 12.00000 |
| EGYPT | 12 | 12.00000 |
| UKRAINE | 11 | 11.00000 |
| DOMINICAN REPUBLIC | 8 | 4.00000 |
| CROATIA | 5 | 2.50000 |
| HONDURAS | 3 | 3.00000 |
| MYANMAR (BURMA) | 3 | 3.00000 |
| SERBIA | 3 | 1.50000 |
| BELGIUM | 2 | 2.00000 |
| CYPRUS | 2 | 2.00000 |
| DJIBOUTI | 2 | 2.00000 |
| AUSTRIA | 1 | 1.00000 |
| BOSNIA-HERZEGOVINA | 1 | 1.00000 |
| NETHERLANDS | 1 | 1.00000 |
| SOUTH AFRICA | 1 | 1.00000 |
| THAILAND | 1 | 1.00000 |
mapCountryData(joinCountryData2Map(disaster_country, nameJoinColumn = "country_code"),
nameColumnToPlot = "death_count",
mapTitle="Earthquakes Total Death")70 codes from your data successfully matched countries in the map
1 codes from your data failed to match with a country code in the map
175 codes from the map weren't represented in your data
mapCountryData(joinCountryData2Map(disaster_country, nameJoinColumn = "country_code"),
nameColumnToPlot = "death_mean",
mapTitle="Earthquakes Death Rate")70 codes from your data successfully matched countries in the map
1 codes from your data failed to match with a country code in the map
175 codes from the map weren't represented in your data
۶. با استفاده از داده لرزه های بزرگ و به وسیله طول، عرض، شدت، عمق مدلی برای پیش بینی تعداد کشته های زلزله بیابید.
مدل glm را باتوجه به نوع دادهها با خانوادهی گاما ایجادمیکنیم و آن را ارزیابی میکنیم و درمییابیم که به عمق و شدت آن بستگی دارد و مکان آن اهمیتی در مدل ندارد پس مدل جدید را براین مبنا مدلسازی میکنیم و ارزیابی میکنیم که نتیجهی قابل قبولی به دست میدهد.
eq_glm = glm(t_death ~ mag + dep + lat + lon, family = Gamma(link = "inverse"), data = disaster)
summary(eq_glm)
Call:
glm(formula = t_death ~ mag + dep + lat + lon, family = Gamma(link = "inverse"),
data = disaster)
Deviance Residuals:
Min 1Q Median 3Q Max
-4.055 -3.273 -2.611 -1.587 7.681
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.158e-03 2.561e-04 4.524 7.90e-06 ***
mag -1.070e-04 2.406e-05 -4.445 1.12e-05 ***
dep 7.511e-06 3.933e-06 1.910 0.0568 .
lat -9.875e-07 2.824e-06 -0.350 0.7267
lon -2.387e-07 4.959e-07 -0.481 0.6304
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 14.9034)
Null deviance: 3965.6 on 429 degrees of freedom
Residual deviance: 3398.8 on 425 degrees of freedom
AIC: 5834
Number of Fisher Scoring iterations: 10
eq_glm = glm(t_death ~ mag + dep, family = Gamma(link = "inverse"), data = disaster)
summary(eq_glm)
Call:
glm(formula = t_death ~ mag + dep, family = Gamma(link = "inverse"),
data = disaster)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.895 -3.270 -2.619 -1.607 7.898
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.223e-03 2.326e-04 5.259 2.30e-07 ***
mag -1.153e-04 2.197e-05 -5.247 2.44e-07 ***
dep 6.468e-06 1.967e-06 3.287 0.0011 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 15.12215)
Null deviance: 3965.6 on 429 degrees of freedom
Residual deviance: 3409.4 on 427 degrees of freedom
AIC: 5832.1
Number of Fisher Scoring iterations: 10
۷. با استفاده از داده worldwide.csv به چند سوال زیر پاسخ دهید. تحقیق کنید آیا می توان از پیش لرزه، زلزله اصلی را پیش بینی کرد؟
از روی دادهها پیشلرزهها را استخراجکرده و روی آنها تعداد و میانگین آنها و قدرت زلزلهی اصلی را درکنار هم جمعکرده و مدلی را استخراج میکنیم که این پیشبینی را صورت دهد که خروجی آن نتیجهی قابل قبولی را گزارش میکند.
worldwide = read_csv("data/worldwide.csv")
worldwide = worldwide %>%
mutate(year = year(time),
month = month(time),
day = day(time),
location = as.data.frame(str_split_fixed(place, ',',2))$V2,
location = ifelse(location == '' , as.character(place), as.character(location))) %>%
select(time, year, month, day, latitude, longitude, depth, mag, location) %>%
group_by(year, month, location) %>%
arrange(time) %>%
mutate(mag_max = max(mag))
worldwide$groupid = worldwide %>% group_indices(year, month, location)
worldwide_temp = worldwide %>%
filter(mag == mag_max) %>%
ungroup() %>%
select(groupid,time_max = time)
worldwide = merge(worldwide, worldwide_temp)
worldwide = worldwide %>%
mutate(pre_past = ifelse(time > time_max, 'pre',ifelse(time == time_max, 'eq', 'past')))
worldwide_pred = worldwide %>%
filter(pre_past == 'pre' | pre_past == 'eq') %>%
ungroup() %>%
select(groupid, mag, pre_past) %>%
group_by(groupid,pre_past) %>%
summarise(avg_mag = mean(mag), count = n()) %>%
ungroup()
wwpm = worldwide_pred %>%
filter(pre_past == 'pre') %>%
select(groupid, count, avg_mag)
wwpmm = worldwide_pred %>%
filter(pre_past == 'eq') %>%
select(groupid, main = avg_mag)
x = merge(wwpm,wwpmm) %>%
select(-groupid)
glm_eq_pred = glm(count~., data = x)
summary(glm_eq_pred)
Call:
glm(formula = count ~ ., data = x)
Deviance Residuals:
Min 1Q Median 3Q Max
-92.62 -15.38 -1.20 10.15 1100.56
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.463 6.533 2.214 0.0269 *
avg_mag -44.653 2.148 -20.790 <2e-16 ***
main 37.713 1.681 22.429 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 2502.604)
Null deviance: 6924258 on 2218 degrees of freedom
Residual deviance: 5545771 on 2216 degrees of freedom
AIC: 23666
Number of Fisher Scoring iterations: 2
۸. گزاره " آیا شدت زلزله به عمق آن بستگی دارد" را تحقیق کنید؟ (طبیعتا از آزمون فرض باید استفاده کنید.)
cor.test(worldwide$depth, worldwide$mag, method = 'spearman')
Spearman's rank correlation rho
data: worldwide$depth and worldwide$mag
S = 4.176e+13, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.2376803
mat_eq = as.matrix(worldwide %>%
select(depth, mag) %>%
mutate(depth = depth + abs(min(depth)))
)
chisq.test(mat_eq)
Pearson's Chi-squared test
data: mat_eq
X-squared = 489570, df = 69011, p-value < 2.2e-16
تست وابستگی بیان میکند که غیرمرتبط بودن این دو متغیر از یکدیگر رد میشود پس میتوان نتیجهگرفت که بین آنها ارتباط وجود دارد و تست استقلال نیز بیان میکند که این دو متغیر نیز مستقل بودنشان از هم رد میشود و وابستگی آنها نیز تایید میشود. پس فرض مستقل بودن این دو متغیر از هم رد میشود.
۹. میانگین سالانه زلزله ها را بر حسب کشور به دست آورید. آیا میتوان دلیلی در تایید یا رد تئوری هارپ ارائه کرد.
long_lat = worldwide %>%
select(longitude, latitude)
long_lat$country = coords2country(long_lat)
worldwide$country_code = long_lat$country
country_code = worldwide %>%
ungroup() %>%
group_by(location) %>%
select(location, country_code) %>%
na.omit() %>%
distinct(country_code, location) %>%
ungroup() %>%
select(location, c_code = country_code) %>%
group_by(location) %>%
slice(1) %>%
ungroup()
worldwide_fixed_country_code = merge(worldwide, country_code)
worldwide_fixed_country_code = worldwide_fixed_country_code %>%
mutate(c_code = ifelse(!is.na(country_code),
as.character(country_code),
as.character(c_code)
),
country_code = c_code
) %>%
select(-c_code)
worldwide = worldwide_fixed_country_code
worldwide_accu = worldwide %>%
group_by(year, country_code) %>%
summarise(eq_count = n(),
mean_mag = mean(mag),
mean_dep = mean(depth))
eq_2015 = worldwide_accu %>%
filter(year == 2015) %>%
select(eq_count)
eq_2016 = worldwide_accu %>%
filter(year == 2016) %>%
select(eq_count)
eq_2017 = worldwide_accu %>%
filter(year == 2017) %>%
select(eq_count)
eq_2018 = worldwide_accu %>%
filter(year == 2018) %>%
select(eq_count)
t.test(eq_2015$eq_count,eq_2016$eq_count)
Welch Two Sample t-test
data: eq_2015$eq_count and eq_2016$eq_count
t = -2.462, df = 109.01, p-value = 0.01538
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-381.34944 -41.19174
sample estimates:
mean of x mean of y
32.98413 244.25472
t.test(eq_2015$eq_count,eq_2017$eq_count)
Welch Two Sample t-test
data: eq_2015$eq_count and eq_2017$eq_count
t = -2.1236, df = 108.3, p-value = 0.03598
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-387.70921 -13.36027
sample estimates:
mean of x mean of y
32.98413 233.51887
t.test(eq_2015$eq_count,eq_2018$eq_count)
Welch Two Sample t-test
data: eq_2015$eq_count and eq_2018$eq_count
t = -1.2653, df = 96.949, p-value = 0.2088
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-214.29194 47.43411
sample estimates:
mean of x mean of y
32.98413 116.41304
با توجه به اینکه پروژهی هارپ در سال ۲۰۱۵ متوقف شد زلزلههای این سال را با سالهای بعدی آن مقایسهکردم و تست t را بر روی آنها اجرا کردم تا تفاوت این سالها را باهم مقایسه کنم. نتایج این تست نشانداد که این سالها بایکدیگر تفاوت معناداری ندارند و این تئوری را میتوان رد کرد.
۱۰. سه حقیقت جالب در مورد زلزله بیابید.
۱۰.۱. نمودار جعبهای زلزلهها در سالها و جایگاه ایران در آن
ww_mag_year = worldwide %>%
group_by(year, country_code) %>%
summarise(mean_mag = mean(mag)) %>%
select(year, mean_mag) %>%
ungroup()
iran_w = worldwide %>%
group_by(year, country_code) %>%
summarise(mean_mag = mean(mag)) %>%
filter(country_code == 'IRN') %>%
select(year, mean_mag) %>%
ungroup()
ggplot(ww_mag_year,aes(x = year, y = mean_mag, group = year)) +
geom_boxplot() +
geom_line(data = iran_w,aes(x = year, y = mean_mag, group = 1, color = 'Iran'), size = 1.5) +
ggtitle("Average Earthquake Magnitude In Worldwide") +
ylab("Magnitude Mean") +
xlab("Year")۱۰.۲. ده زلزلهخیزترین کشورها
country_reg = country_code %>%
group_by(c_code) %>%
summarise(reg_num = n()) %>%
ungroup() %>%
select(country_code = c_code, reg_num)
ww_mag_count = worldwide %>%
group_by(country_code) %>%
summarise(mean_mag = mean(mag), count = n())
ww_mag_count = merge(ww_mag_count, country_reg)
top_10 = ww_mag_count %>%
mutate(score = (count*mean_mag)/reg_num) %>%
arrange(-score) %>%
slice(1:10)
knitr::kable(top_10, "html")| country_code | mean_mag | count | reg_num | score |
|---|---|---|---|---|
| JPN | 4.531295 | 2486 | 1 | 11264.80 |
| PNG | 4.576130 | 2191 | 1 | 10026.30 |
| NZL | 4.537366 | 2058 | 1 | 9337.90 |
| PRI | 2.886450 | 2972 | 1 | 8578.53 |
| PHL | 4.531641 | 1438 | 1 | 6516.50 |
| IDN | 4.462053 | 4277 | 3 | 6361.40 |
| TON | 4.628141 | 1329 | 1 | 6150.80 |
| VGB | 2.986832 | 1632 | 1 | 4874.51 |
| VUT | 4.580264 | 1059 | 1 | 4850.50 |
| RUS | 4.425916 | 1092 | 1 | 4833.10 |
۱۰.۳. نحسی عدد ۱۳ در وقوع زلزله در روز ۱۳ام.
thirteen = worldwide %>%
filter(day == 13)
not_thirteen = worldwide %>%
filter(day != 13)
t.test(thirteen$mag, not_thirteen$mag, alternative = "greater")
Welch Two Sample t-test
data: thirteen$mag and not_thirteen$mag
t = 6.5195, df = 2091, p-value = 4.407e-11
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
0.0961231 Inf
sample estimates:
mean of x mean of y
3.873674 3.745097
این تست رد میشود و نمیتوان نتیجهگرفت که روز ۱۳ام روز نحسی برای زلزله است.